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ABSTRACT 

This paper presents an adaptive method that can he used to improve the finite element solution of a linear or 
non-linear PDE by indicating the local as well as global errors. This wavelet-based adaptation is based on hierarchical 
finite element error estimation scheme. The finite element solution of engineering problems is transformed into the 
multilevel decomposition of wavelet space. The functions of wavelets coefficient are used as an a posteriori error 
indicator. The error estimations in the present setting do not have any problem due to the complex domain of engineering 
structures and is applicable in some small region as well as the complete domain. The proposed technique is quite 
efficient for linear problems like structures with discontinuity as well as non-linear ones such as contact problem. 
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1. INTRODUCTION 

There is no analytic theory for mesh refinement for most of the boundary value problems (BVP). 
Therefore, the adaptive mesh cannot be generated to improve the solution on the basis of any prior information. 
A posteriori error estimate is a key ingredient for improving the accuracy of the solution in the successive iteration. 
The adaptive mesh generation process depends on the successful application of a posteriori error estimation 
techniques which is not only a tool for grid optimization but also provide the estimate of the error in the current 
numerical solution. There exists a vast literature on various error estimation techniques [1-12] which are recently 
developed. Most of the estimators use local quantities such as jump of derivatives across the interface between two 
adjacent elements[13-18]. The residuals due to local quantities are used for estimating upper and lower bounds of 
the global error in terms of energy norm. The quality of estimated error by using canonical error-estimation 
techniques depends on the initial mesh which is often overestimated. 

The hierarchical approach, which is a multi-scale expansion[19-20]of corresponding FEM, is often used 
for a posteriori error estimate due to its simplicity and effectiveness[21-22]. The hierarchical basis functions span 
the space of traditional basis functions of finite element method, and it has desirable properties for fast iterative 
solvers. A posteriori error in the hierarchical approach is represented as a convergent sequence of errors of the 
Hilbert space. In other words, unlike the other classical error estimators which use total error within an element as 
an indicator, the higher level component is used as an indicator in this scheme. It is noted that such a posteriori 
error estimators are useful in identifying local as well as global error. 
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The hierarchical basis functions are defined as “lazy wavelets”[23]. The hierarchical scheme has been applied to 
various types of partial differential equations, but wavelet-based techniques are not used frequently. Classical wavelets 
have many desirable properties such as orthogonality and vanishing moments which do not exist in a hierarchical basis. 
Bertoluzza used wavelet-based error estimator for wavelet Galerkin method. A good mathematical theory for error 
estimation and adaptive solution of PDEs using wavelet basis are developed[25,26]. Implementation of these 
wavelet-based methods is difficult on the wide range of complex boundary value engineering problems. Sudarshan[27] et 
al have used wavelet-based a posteriori error estimation for engineering structures, but this method requires customized 
wavelets. It is very difficult to develop customized wavelets, which should be as efficient as FEM for various types of 
partial differential equations. 

The present work is based on the hierarchical error estimation technique proposed by Bank[28]. It combines the 
strength of a wavelet transform and hierarchical error estimation technique and finds the error in the finite element 
solution. In this method, the solution of the finite element method is transformed into the multi-resolution wavelet space 
using a linear transformation. For a single variable system, the detail coefficients are used as a local error indicator but 
some functional of detail coefficients is to be used for multiple variable problems. In the case of variational inequalities, we 
solve a problem that satisfies the assumptions made for the linear case and calculated the error bounds defined by 
Kornhuber et al[29]. 

2. WAVELET-BASEDERROR ESTIMATION 

Consider a partial differential equation Lu = / in Q and boundary conditions Bu — g . The problem is to find 
the exact solution U G TL . The variational form of the equation can be expressed as 

a(u,v) = f(v) (1) 

Forall VE TL , where TL is an appropriate Solve space,«(.,.) is a positive definite bilinear form, and / (.) is linear 
functional. The energy norm associated with a( .,.) can be defined as: 

II i|2 / x 

p =a(u,u ). (2) 

The best available option to solve (1) is the finite element method, and the current solution is assumed as the 
highest resolution for the wavelet space Uj eVj CzTL . This approximate solution process introduces numerical error 

e = 11 - llj . (3) 

Obtaining the exact error £is as difficult as finding the exact solution. Thus, various indicators for a posteriori 
error estimations are developed. If the estimated error is larger than a predefined tolerance, then a re-meshing is used on 
the finite element grid either selectively on element with large errors or on the whole domain (uniform refinement). 
This is essential to form a strategy of adaptive mesh for the improvement of the solution in the successive iteration. 

A brief introduction of wavelets will be useful for this error estimation technique. The concept of multi-scale 
analysis is to interpolate an unknown field at a coarse level with the help of so-called scaling functions. Any improvement 
to the initial approximation is achieved by adding ‘details’ provided by new functions known as wavelets. 
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A multiscale analysis uses multire solution properties of wavelets[30]. Each subspace V■ is spanned by a set of scaling 
functionj^. k (x),Vfce Zj . The complement of Vj in Vy +1 is defined as subspace Wj such that \L +1 = Vj ©Wyfor all 
j E Z . Space Vj +1 can be decomposed in a consecutive manner as: 

V j +i=Vo @ Wo ® Wi ® Wr ''® Wj . (4) 

The basic functions inWjare called wavelet functions and are denoted byl//j k . LetV ; d L 2 (£l)be a finite 
dimensional subspaces, and consider the approximate solution U ■ e Vj such that fljw^vj = f (v) for allveVy. Its 


solution holds the approximation property 


U-Uj 


— inf ||l/ — v||. Similarly, for Vj d Vj +l d L 2 (fl) Galerkin solution 


Uj +l e V 7+1 will satisfy ^|w ;+1 ,vj = /(v)for allve \L +1 . Here U—Uj —>0 , as j —> oo . 


Let us assume that as compared W • the next scale refined solution U . +] is closer to U . Thisis definedin terms of the 


saturation assumption 


u-u j+l 


U-Uj 


(5) 


Where < 1. The characteristic length of a mesh element h —> 0 as j —> °° . Lor a higher degree of 
approximation in the space Vy +1 , Bank [21]estimated dj = 0{h r j for someT>0. Lor V^eV^ the unique 
decomposition V . +1 = Vj + Wy, where V ; - E Vy and Wj E V/ , we assume a strengthened Cauchy inequality for the 
decomposition, 

( 6 ) 


lalv^W: ) 

<X 

V; 

W: 

w’ j) 


J 

J 


Where X < 1 is independent of j. 

Instead of (3), Bank [21]used two-scale error in the multi-resolution approach 


r j= U j +l -Uj. 


(7) 


Where r • E W■ , it means that two-level error is the projection of true error onto the wavelet space. It is proved 


[21]that 


c x U~Uj < Vj <c 2 U~Uj . 


( 8 ) 


Where C x , C 2 >0 are constants and independent of the multi-resolution level j. The orthogonality property [26] 

2 



2 


2 


U-Uj 

— 

u-u j+ 1 

+ 

Uj+l - Uj 


(9) 
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The saturation property is used to get the lower bounds C x — |l — j. The two scale error will always be less than 

. Therefore, the two-level error is a good estimate of the actual error. Since w. is the 


the total error, i. e. 

r- 

< 

u-u • 


j 


J 


projection of r hence 


W 


7+1 


< 


r ;+i 


< 




< 


< 


U-Uj 


( 10 ) 


The two consecutive levels of FEM refinements can be used to calculate two scale error r . +] or T- . 

This method is expensive to estimate error. Therefore FE solution is transformed using discrete wavelets. The norm of 
some functional calculated using wavelet coefficients shows a sudden jump in some areas which indicate the need for mesh 


refinement in that zone. The cost of computing 




or 


w j 


is considerably less than the two scale error. 


Theorem 

Let V j+1 = Vj © Wj , strengthen Cauchy inequality and saturation property holds then 


(l-<^ 2 )(l-X 2 )|| u—Uj 


< 


w j 


< 


U-Uj 


Proof 


The proposed posteriori error estimate Wy E Wj is the error of the solution Uj E V- such that 
a(w ; -,v) = /(v)-a(M 7 -,v) V veWj 
ox,a[u-Uj -Wj,v) = 0 V veWj. 


( 11 ) 


A /V 

Substitution of V = Wj in the preceding equation gives right inequality of the theorem. Let Uj +1 = Uj + Wj where 


Uj E Vj and Wj E Wj then 


> 


> 


j + l~ u j 

— 

u j- 

u- —u- 

2 

+ 

W: 

J J 


J 

u —ll , 

(x 

W: 

J J 

\ 

J 


- 2 % 


u i~ u j 


w j 


- 2 % 


W. 


+ 


w, 


>(1-X 2 )| 


w. 


( 12 ) 


Impact Factor (JCC): 7.6197 


SCOPUS Indexed Journal 


NAAS Rating: 3.11 


































































































A Customized Wavelet-Based Hierarchical A Posteriori 
Error Indicator for Linear and Nonlinear PDES 


321 


a 


By using orthogonality relations a | U — Uj ,Vj — 0 V VE V j and a[u — Uj + i ,VJ — 0 V VE Vj + \ , we get 
| u j + 1 — , V J = 0 V veVj . Substitution of V = — U -, gives 




= a ( Uj 


U„W: . 


7+1 “ 7 ’ 7 


(13) 


On subtracting (11) from a[u— My + j,vj — 0 V V’G Vy + |, we get| 
With V = Wj , we get 

a(u j+l -Uj,Wj) = a(wj,Wj). 

Using (13) and (14), we have 

\\u j+l - Uj f = a (Wj,Wj). 

Using (6) in the preceding equation, it can be shown that 


U j+X -Uj-Wj 



2 



u j+l~ u j 

< 

w j 

w j 


From (12) and (15) we obtain 

*M 2 )N- 


w 


On substituting (15) in (9), we get 



2 


2 



U~Uj 

— 

+ 

S 

1 

a 

+ 

W 7 



Using saturation assumption, we obtain 


U-U • 


<^ 2 


U—Uj 


+ 


W. 


W. 


Combining this with (16), we have 


(l-^ 2 )|| U-Uj 


< 


1-X 2 


w j 


.,v) =0 V veWj. 


(14) 


(15) 


(16) 


3. ESTIMATION OF WAVELET COEFFICIENTS 

Now we have to determine the error in a computationally efficient way. In this paper, we used the functional of 
wavelet coefficients. A function M- +1 E L}(R) is approximated by its projection P^Uj onto the space Vj and the 

projection on W- as Q J u ■, 
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P J+ 'u , +l = P j Ui +Q j u,. 

7 +1 J ^ J 

The coefficient vector Uj = {Uj q,. ,Uj } T of some scaling function basis forms the projection P j Uj and 

Wj = {Wj 0 ,., W. m( n } T of some wavelet forms the matrixQ^.. The approximate solution /i.gV. czV can be 


jMj) 


represented in a multi-resolution format consisting of scaling functions and compactly supported wavelets up to a 
resolution of j 


7-1 


U M = Z U O,k0O,k + Z Z r j,mVj,, 

k j =0 m 


(17) 


Here k and r . m are a projection of function U on scaling function (j) and wavelets \j/ respectively. We take 
to be Lagrange finite element functions and Ijf as hierarchical finite element bases [31]. 

Considering only one level we have U. 


7+1 


7+1 


: E«iA-+Zw, 


(18) 


The error W. becomes(from equation (7)) 

II r. 11=11 Y u iAt+ll r j,sy / i,s-Y u jAt 

t S t 

="L o.y'j," 


(19) 


To calculate the detail coefficient, we use the multi-resolution format of the trial function in the Galerkin 
formulation 


j -1 


«(Z U oA ,k + Z Z ’ V h) = l ( V h) 

k j =0 m 

We get a multi-scale stiffness matrix in the form 


( 20 ) 


cs C/ 

J o,i 

O/ 'i’ nD 

Vi u i 

c/ * 7 c/ . 7 

■ 0j —1 ^1,7-1 


S 1 

^Oj-l 

c/ . 

°lj-l 

nP 

u j-i 


- 


r s 

1 M o 


f ) 

r o 

= 

f D 

7 0 

1 - 

Oi 

i_ 


1 

1_ 


( 21 ) 


Here the sub-matrices V and D.° are defined as 
S s = a{(f> 0 ,(f> Q ) where (f> 0 = [<t> 0k } 


( 22 ) 
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Dj = a(y/ r y/.) where (/( ={ y/ i k } 

(23) 


They contain interactions of scaling -scaling and wavelet-wavelet functions spanning sub-spaces of the same 
scale. Sub-matrices S 1 are defined as 

7-1 


Sij-i= a Wi- (24) 

They contain connection coefficients of interactions between scaling and wavelet functions. Vectors f u and w 
represent the decomposition of load vector and the unknown solution over different scales of resolution. By solving the sub 

system D d r t = f t D we can calculate the two-level error/details at a certain level/. But as the multi-resolution matrix is 

coupled over scales, we need to eliminate the interaction sub-matrices. The interaction sub-matrices can be eliminated by 
customizing the wavelet basis \j/ [27]such that 

S[j_ x = 0 (25) 

So, the system becomes decoupled over scales. Scale decoupling of the multilevel system splits it into sub- 
systems that can be solved individually. Such wavelets that can decouple the multi-scale matrix are created using the 
second generation wavelet techniques, i. e. lifting scheme and stable completion [27]. 

Vi,j ~ ^ "*" ^ S i,m,jfii,m @ 6 ) 

k m 

Here Xj/; y an d 0,- y are the j th wavelet and the Lagrangian basis function, respectively at the i th level of resolution 
g and S coefficients are obtained by solving the system 

K .X , = 0 (27) 

i,j i,j v 7 


In this equation, K t y is a block interactions sub-matrix with respect to the weak form of the desired PDE operator 
and is defined as 

K i,}' = ’ Z 8i,k,f0i+l,k + Z S i,m,jAm ) ( 2 ») 

k m 

This sub-matrix occurs naturally in the multi-scale stiffness matrix of equation(21), and is composed of the inner 
product between k + m Lagrangian basis functions, the wavelet function y/. . and all Lagrangian basis functions that have 

overlapping local support on the j ' th level. X . is the matrix with n columns of k+m coefficients g and S , and n is the 
number of resulting wavelets for the level j '. 


The algebraic form of the linear variational problema|^y ,vj = / (v) is 


Find Uj e Vj : KjUj = f j (29) 

Using the multi-resolution structure of trial function, we get a multi-level stiffness matrix (equation 21). Using the 
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combination of equations (20),(27),(28) some multi-scale system matrix (equation 21) can be block diagonalized 



1 

^ o 

3 

1 _ 


i 

i_ 

. . 

r 

'o 

= 

f D 

j 0 

. . D%_ 

i - 

i_ 


f D 


The first sub-system 

« = fo 


(30) 


(31) 


will produce the initial approximate solution of low accuracy, and is same as the finite element solution obtained 
in a coarsely discretized mesh. Solution of subsystem 

= / 0 D (32) 


results in details belonging to wavelets of scale 1. The details are then combined with the coarse scale solution, 
providing the improved solution. We continue the refinement process until the desired accuracy is reached. 

The detail coefficient computed at each level of resolution is the difference between two successive solutions, so 
it iscalled as a two-level error. Use of wavelet-based error estimation doesn't require post-processing of results. The other 
methods of error estimation require an evaluation of specific equations on an element by element basis. 

4. NUMERICAL EXPERIMENTS 


We developed the wavelet-basedcost-effective error indicator. To test the performance of proposed error 
estimators, we compare it with exact error and two scale error indicator for standard problems. In the first example, we 
consider the linear initial-boundary value hyperbolic problem: 

^L-y^L+x^- = 0, -l.2<x, y<\2, t> 0, 
ot ox oy 


and initial conditions u(x, y,0) = 


J0 if (x-0.5) 2 +1.5j 2 >0.0625 

|l-16((x-0.5) 2 +1.5y 2 ) otherwise. 


The Dirichlet boundary conditions on all four sides u(l.2,y,t) = u(—1.2,y,t) = u(x,—1.2,t) = u(x,1.2,t) = 0 are 
considered. The solution to this problem is a moving elliptical cone that rotates counterclockwise direction around the 

where 


origin with a period 71 . It can be written in the form 

C = l-16[(xcost + ysint-0.5) 2 +1.5(ycos^ + xsin0 2 ]. 


u(x,y,t) - 


if C<0 
if C > 0, 


The equation is solved by using finite element method for various level of uniform mesh. AtO < t < 3.2, the 
New mark time integration scheme is used with a time step A t = 0.00625 s. 
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Level of refinement 

Figure 1: Comparison of Exact, Two Scales, and Proposed Error Estimators 


Figure 1 shows the errors 


U—Uj 


r, 



j 


and the proposed estimator 


w ,- 


with the increasing degrees of freedom. 


Initially at j = 1, a uniform grid of 20 X 20 is used. The grid is doubled at every successive level j for multi-scale 
wavelets. The figure shows the convergence of exact error as the degrees of freedom increases. It can be observed that the 
two-level indicator and the proposed indicator show the similar trend of convergence of global error with the level of 
refinement. This ensures the reliability of the proposed wavelet-based estimator. 



Figure 2: One-Quarter of the Discretized Rectangular Plate with a Hole 
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Figure 3: Strain Energy using Two Scale Difference 

To demonstrate the effectiveness of the two-level error indicator and proposed error indicator for local error 
estimation, both the methods are tested on a standard 2D linear elasticity problem. The zone of mesh refinement for a plate 
with a hole under longitudinal tensile loading is very well known. Both the error indicators are tested by considering one- 
fourth of the plate as shown in Figure 2. The regions of interest in these problems are the elements in the neighborhood of 
the circular hole where stress concentration is very high. The functional used for two-level error indicator is the difference 
of strain energy at two consecutive levels of finite element solutions. The peaks in the Figure 3, indicates the elements with 
high local error i. e., these elements are to be refined. For the proposed error indicator, the strain energy is calculated by 
using wavelet coefficients and shown in Figure 4. As in the previous example, we see the similar peaks at high-stress 
concentration elements, but the proposed error indicator is much more economical than the two-level error indicator. 
Convergence of stress is used to stop the refinement process. 



Figure 4: Strain Energy using Wavelet Coefficient 
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As an on-linear example, we consider the contact problem that is a variational inequality of the first kind. The 
governing differential equation of contact can be written in the strong form as 

-V (C : £(u )) = / in Q ( 33 ) 

u = 0 on T d , 
crn = t on V N , 


u n < g andn^cm > 0 on T c . 

Here T represents the boundary of the domain, with F c , F d , F n representing the contact, Dirichlet, and Neumann 

boundaries respectively( F c V d V n = T) with g defined as gap function and nis the unit normal along the contact 

boundary. 

The inequality is written in the bilinear form as 

a(u,v- u) > Z(v- u) Vve K (34) 

K is defined as a closed and convex set of permissible functions K = {y E fly (^j )l v.n < g) 

We discretize (34) using </>. e V such that a finite-dimensional dimensional problem is obtained 


" £K i --a(Uj,Vj- u J )>I(v i -Uj) Vvj e K. 


(35) 


With K ^ given as 


K / = v e n j Vj(Qy) / r v.tijij ds < j r (g - Uj. n) Hjds, V// e M } 


Here M is the cone of Lagrange finite element function [32] 

The two-level error satisfies 

a(r. , v - r .) > p(v - r .) Vv e K j+l , (36) 

and residual p is defined as 


p = l(v)-a(u.,v.) (37) 

With the help of triangle inequality and saturation assumption 

(l-X) II u - u. Il<ll r. Il< (1 + X) II u - u. II 

The contact problem is a nonlinear problem so it is solved iteratively. To calculate the two-level error we compute 
the stiffness matrices of the contact problem at a load step i . For the variational inequality, the linearized algebraic form of 
equation (36) is written as 
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Find Uj E K .: (K'. + K' mj )Au'. 


f ; 


(38) 


Where, K . is the consistent global matrix of dimension A f 1 . X A f 1 . , and K l m ■ is the geometric tangential stiffness 
matrix which is a diagonal matrix of dimension X dim( JVj ) . A u\ is the incremental displacement at load step i. Here 

J\f\ is the total number of mesh nodes which are not on contact boundary A f C is the number of nodes on the contact 
boundary and AA is the total number of nodes at mesh resolution j. Using the multi-resolution and second generation 
wavelet customization, the algebraic equation at each iteration step has the form 


A 


Ai<q 


1 

O 

Ci 

_1 


Au s 0 


~/o " 

D? . . 


Ar o 

+ 




= 

/o° 

■ ■ D U. 

i 

Ah. 

i 

• <+_ 

i 

Ar. , 

L 7-1 J 

i 

_A_ 


We solve a small sub-system to calculate the wavelet detail 


(39) 


(40) 


A Hertz contact problem with the known analytical solution is considered to check the efficiency of our error 
estimator for variational inequalities[33]. The contact problem involves a frictionless contact condition between two 
infinitely long cylinders. The material is linearly elastic and isotropic, and deformation is considered infinitesimal with 
plane strain condition. For contact condition between a cylinder and a plane surface shown in figure 5, we assumed large 

values of radius R 2 and Young’s modulus E 2 for the rigid surface. The Hertz analytical solution for the normal contact 
stress (jH and the length b are given as [33] 




7tb 


l 2R?q(l-V?) 
E x n 


b-2 


We assumed Poisson’s ratio V x = 0.3, the radius of the cylinder R { = 8mm, Young’s modulus E } = 80 GPa and the 

force applied to the structure F = 100 N. We compute the distributed load q = 625 N/m and the contact region b = 0.68077. 
The plane surface is taken as a rigid surface with a constrained lower edge. For the initial contact, closed gap condition is 
taken at the symmetry point of the cylinder. 

Figure 5 shows the computational domain and a reduced model for simulation. The number of total nodes in the 
initial mesh is 150. Linear Dirichlet boundary condition is considered for non-contact points. The wavelet coefficients are 
estimated at each iteration step for local refinement. The normalized detail coefficients are measured in the energy 
norm[31]. 
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Figure 5: Two Cylinders in Contact and the Rigid-Deformable Approximation 



i 


Figure 6: Amplitude of Normalized Wavelet Coefficients 

The quantity of interest in this problem is maximum contact stress. The zone of interest in this problem is the 
elements around the neighborhood of contact point as the point of maximum stress lies on the contact boundary. The 
wavelet-based adaptive scheme is compared with uniform /z-refinement of the computational mesh. An idealized 
equilibrated force system is applied at the contact interface. Figure 6 shows the distribution of details at refinement level 3; 
the wavelet coefficients show a peak over the elements with a maximum local error. The mesh is refined at the grid 
showing the high amplitude of wavelet detail. 


x 

Figure 7: Von Mises Stress Distribution 

Figure 7 shows the von Mises contact stress distribution at refinement level 3, and it can be observed the 
maximum contact stress is coming at the point of contact. Figure 8 shows a comparison between wavelet-based adaptive 
refinement and uniform refinement. it can be obsereved that for the same level of error reduction wavelet-based error 
indicator is more economical than uniform refinement. 
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Figure 8: Comparison of Uniform vs. Wavelet-based Refinement 

5. CONCLUSIONS 

Accuracy, reliability, and computational cost are crucial issues of any efficient error estimator to design adaptive 
algorithms for the numerical solution of partial differential equations. The wavelet-based a posteriori error indicator is 
proposed to identify the location of grid refinement at a low computational cost to achieve an error bound solution for 
engineering problems. In this approach, the wavelet transform is used for error indicators that are easier and much simpler 
to implement than those available in the literature. Furthermore, it is locally efficient and can be used in complex 
engineering analysis with any finite element subspaces. Comparisons of the proposed estimator with exact error and two 
scale error with the help of numerical experiments show that it is as good as other error estimators. The key feature of the 
method is an automatic algorithm for any finite element grid structure. In situations where it might be impossible to 
directly verify the hypothesis, the findings of this study facilitate future researchers to advance the error estimation 
techniques. There are several areas in wavelet-based error estimator still needs research. This includes the selection of the 
most suitable wavelets for complex geometries with the large class of practical applications. 
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